Coweeta Hydrologic Laboratory Watershed 32 using Static, SSURGO, and RSS soil inputs.

The following R Markdown script prepares and runs the Regional Hydro-Ecologic Simulation System (RHESSYS) within the R environment. Modified from RHESSysIOinR example scripts created by Will Burke and Ryan Bart. Includes CLHS code snippets and SDA code from Dylan Beaudette.

Windows Subsystem via Linux must be installed before RHESSys will run in a Windows environment. GRASS8 must be installed to run the spatial preprocessing portion of the script.

Code has been modified to run RHESSys 7.4 and 5.18. 7.4 vegetation processing differs slightly from 5.18. Various template differences exist between the two versions of RHESSys and template files must be carefully reviewed if moving between versions. This script is currently set up to run RHESSys 7.4

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
#install.packages("devtools")
#devtools::install_github("ropensci/FedData")

library(aqp)
## This is aqp 1.41
library(daymetr)
library(clhs)
## Registered S3 methods overwritten by 'tibble':
##   method     from  
##   format.tbl pillar
##   print.tbl  pillar
library(EcoHydRology)
## Loading required package: operators
## 
## Attaching package: 'operators'
## The following objects are masked from 'package:base':
## 
##     options, strrep
## Loading required package: topmodel
## Loading required package: DEoptim
## Loading required package: parallel
## 
## DEoptim package
## Differential Evolution algorithm in R
## Authors: D. Ardia, K. Mullen, B. Peterson and J. Ulrich
## Loading required package: XML
library(elevatr)
library(ezknitr)
library(FedData)
## 
## Attaching package: 'FedData'
## The following object is masked from 'package:operators':
## 
##     %>%
library(foreach)
library(foreign)
library(ggplot2)
library(ggpubr)
## 
## Attaching package: 'ggpubr'
## 
## The following object is masked from 'package:operators':
## 
##     %>%
library(hydroGOF)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'hydroGOF'
## The following object is masked from 'package:topmodel':
## 
##     NSeff
library(imputeTS)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attaching package: 'imputeTS'
## The following object is masked from 'package:zoo':
## 
##     na.locf
## The following object is masked from 'package:operators':
## 
##     %>%
library(knitr)
library(latticeExtra)
## Loading required package: lattice
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(Metrics)
## 
## Attaching package: 'Metrics'
## The following objects are masked from 'package:hydroGOF':
## 
##     mae, mse, rmse
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following objects are masked from 'package:aqp':
## 
##     metadata, metadata<-
library(rasterVis)
library(readxl)
library(reshape2)
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: C:/Users/Carlos/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Carlos/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(rgrass)
## GRASS GIS interface loaded with GRASS version: GRASS 8.2.1 (2023)
## and location: Coweeta1
library(rmarkdown)
library(Rmisc)
## Loading required package: plyr
## 
## Attaching package: 'plyr'
## The following object is masked from 'package:ggpubr':
## 
##     mutate
library(RHESSysPreprocessing)
library(RHESSysIOinR)
## 
## Attaching package: 'RHESSysIOinR'
## The following object is masked from 'package:RHESSysPreprocessing':
## 
##     read_world
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
## 
## Attaching package: 'sf'
## The following object is masked from 'package:operators':
## 
##     %>%
library(soilDB)
library(sp)
library(stats)
library(stringi)
library(tactile)
library(terra)
## terra 1.5.21
## 
## Attaching package: 'terra'
## The following object is masked from 'package:rgdal':
## 
##     project
## The following object is masked from 'package:knitr':
## 
##     spin
## The following object is masked from 'package:zoo':
## 
##     time<-
## The following object is masked from 'package:ggpubr':
## 
##     rotate
## The following object is masked from 'package:ggplot2':
## 
##     arrow
library(testthat)
## 
## Attaching package: 'testthat'
## The following object is masked from 'package:terra':
## 
##     describe
## The following object is masked from 'package:operators':
## 
##     %>%
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble  3.0.3     v dplyr   1.0.5
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## v purrr   0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x forcats::%>%()           masks stringr::%>%(), dplyr::%>%(), purrr::%>%(), tidyr::%>%(), tibble::%>%(), testthat::%>%(), sf::%>%(), imputeTS::%>%(), ggpubr::%>%(), FedData::%>%(), operators::%>%()
## x purrr::accumulate()      masks foreach::accumulate()
## x dplyr::arrange()         masks plyr::arrange()
## x terra::arrow()           masks ggplot2::arrow()
## x lubridate::as.difftime() masks base::as.difftime()
## x dplyr::combine()         masks aqp::combine()
## x purrr::compact()         masks plyr::compact()
## x dplyr::count()           masks plyr::count()
## x lubridate::date()        masks base::date()
## x readr::edition_get()     masks testthat::edition_get()
## x tidyr::extract()         masks terra::extract(), raster::extract()
## x dplyr::failwith()        masks plyr::failwith()
## x dplyr::filter()          masks stats::filter()
## x dplyr::id()              masks plyr::id()
## x terra::intersect()       masks raster::intersect(), lubridate::intersect(), base::intersect()
## x purrr::is_null()         masks testthat::is_null()
## x dplyr::lag()             masks stats::lag()
## x latticeExtra::layer()    masks ggplot2::layer()
## x readr::local_edition()   masks testthat::local_edition()
## x dplyr::matches()         masks tidyr::matches(), testthat::matches()
## x dplyr::mutate()          masks plyr::mutate(), ggpubr::mutate()
## x dplyr::rename()          masks plyr::rename()
## x dplyr::select()          masks raster::select()
## x lubridate::setdiff()     masks base::setdiff()
## x dplyr::slice()           masks aqp::slice()
## x dplyr::src()             masks terra::src()
## x dplyr::summarise()       masks plyr::summarise()
## x dplyr::summarize()       masks plyr::summarize()
## x terra::union()           masks raster::union(), lubridate::union(), base::union()
## x purrr::when()            masks foreach::when()
library(viridisLite)
library(dplyr)


setwd(system.file("extdata/", package = "RHESSysIOinR"))
 
rh_ver = dir(path = "modelfiles/rhessys/", pattern = "^rhessys\\d+",recursive = F)

getwd()
## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"
rh_path = file.path("modelfiles/rhessys", rh_ver)

## Load in necessary files
filedir <- ("C:/Users/Carlos/Documents/GitHub/RHESSys74-R-Markdown/")

watershedshapefile <- paste0(filedir, "input_files/shapefiles/watersheds/Coweeta_Hydrologic_Laboratory.shp")
streamshapefile<- paste0(filedir,"input_files/shapefiles/streams/StreamsCaldwellClipped.shp")
weirshapefile<- paste0(filedir,"input_files/shapefiles/weirs/Subbasinoutlets.shp")
roadshapefile <- paste0(filedir, "input_files/shapefiles/roads/UpdatedRoads4Reprojected32617.shp")

CWRG12 <- paste0(filedir,"input_files/precipitation/RDS-2017-0031/Data/RG12_daily_1942_2021.csv")
CWRG20 <-paste0(filedir,"input_files/precipitation/RDS-2017-0031/Data/RG20_daily_1962_2021.csv")

CS01 <- ("C:/Users/Carlos/Desktop/ORISE/Climate Data/CW/Air Temperature/CS01 OISHI/cs01_hourly_Tair_2011_2023.csv")
CS01LT<- ("C:/Users/Carlos/Desktop/ORISE/Climate Data/CW/Air Temperature/RDS-2015-0049/Data/cs01_daily_airtemp.csv")
CS21<- paste0(filedir,"input_files/temperature/RDS-2015-0042/Data/cs21_daily.csv")
CS28<- paste0(filedir,"input_files/temperature/RDS-2015-0042/Data/cs28_daily.csv")

newsmlocations<-  paste0(filedir,"input_files/shapefiles/soil_moisture/CoweetaNRCSsmKML.kml")

Plot and Filter Observed Soil Moisture

This data has not been pre-prepared and must be cleaned

“Legacy” Soil moisture datasets

Coweeta Long Term Ecological Research Program and P. Bolstad. 2019. Continuously measured soil moisture, soil temperature, and air temperature from stations in Watershed 32, Coweeta Hydrologic Laboratory ver 19. Environmental Data Initiative. https://doi.org/10.6073/pasta/e47b36d4a18b8009b2a786acebeb9cba (Accessed 2024-06-17).

setwd(system.file("extdata/", package = "RHESSysIOinR"))


smws32_1<- read.delim(paste0(filedir,"input_files/soil_moisture/knb-lter-cwt.1308.19/CWT_132_1_0.txt"), header = T)
smws32_2<- read.delim(paste0(filedir,"input_files/soil_moisture/knb-lter-cwt.1308.19/CWT_232_1_1308.txt"), header = T)
smws32_3<- read.delim(paste0(filedir,"input_files/soil_moisture/knb-lter-cwt.1308.19/CWT_332_1_1308.txt"), header = T)

smws32_1$Date<- as.Date(smws32_1$Date)
smws32_2$Date<- as.Date(smws32_2$Date)
smws32_3$Date<- as.Date(smws32_3$Date)

smws32_1 <- smws32_1[,!sapply(smws32_1,is.character)]
smws32_2 <- smws32_2[,!sapply(smws32_2,is.character)]
smws32_3 <- smws32_3[,!sapply(smws32_3,is.character)]

#range(smws32_3$smois30A)

smws32_1mean<- aggregate(smws32_1, by = list(smws32_1$Date), FUN = function(x) mean(x, na.rm = TRUE))
smws32_2mean<- aggregate(smws32_2, by = list(smws32_2$Date), FUN = function(x) mean(x, na.rm = TRUE))
smws32_3mean<- aggregate(smws32_3, by = list(smws32_3$Date), FUN = function(x) mean(x, na.rm = TRUE))

smws32_1mean <- smws32_1mean[order(smws32_1mean$Date),]
smws32_2mean <- smws32_2mean[order(smws32_2mean$Date),]
smws32_3mean <- smws32_3mean[order(smws32_3mean$Date),]

#Filter Soil Moisture Datasets

smws32_1mean  <- smws32_1mean[smws32_1mean$smois30A>'0',]
smws32_3mean  <- smws32_3mean[smws32_3mean$smois30A<='0.45',]
smws32_3mean  <- smws32_3mean[smws32_3mean$Date<='2018-05-15',]


plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", main = "Soil Moisture 30cm Site 1", xlab = "Date", ylab = "Soil Moisture 30A")

plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", main = "Soil Moisture 30cm Site 2", xlab = "Date", ylab = "Soil Moisture 30A")

plot(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "red", main = "Soil Moisture 30cm Site 3", xlab = "Date", ylab = "Soil Moisture 30A")

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Observed Soil Moisture WS32 30cm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}

smws32_1mean  <- smws32_1mean[smws32_1mean$Date>='2017-01-01',]
smws32_2mean  <- smws32_2mean[smws32_2mean$Date>='2017-01-01',]
smws32_3mean  <- smws32_3mean[smws32_3mean$Date>='2017-01-01',]

smws32_1mean  <- smws32_1mean[smws32_1mean$Date<='2018-05-15',]
smws32_2mean  <- smws32_2mean[smws32_2mean$Date<='2018-05-15',]
smws32_3mean  <- smws32_3mean[smws32_3mean$Date<='2018-05-15',]

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 30cm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}

{plot(density(smws32_1mean$smois30A), ylim = c(0,14), xlim = c(0,0.35), col = 'red', main = "Density of Observed Soil Moisture WS32 30cm")
lines(density(smws32_2mean$smois30A), col = 'blue')
lines(density(smws32_3mean$smois30A), col = 'black')
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 1")
lines(smws32_1mean$Group.1,smws32_1mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_1mean$Group.1,smws32_1mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_1mean$Group.1,smws32_1mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1-30A","SM1-30B","SM1-60A","SM1-60B"), col=c("red","orange","blue","black"), lty=1)}

{plot(density(smws32_1mean$smois30A), ylim = c(0,14), xlim = c(0,0.35), col = 'red', main = "Density of Observed Soil Moisture WS32 Station 1")
lines(density(smws32_1mean$smois30B), col = 'orange')
lines(density(smws32_1mean$smois60A), col = 'blue')
lines(density(smws32_1mean$smois60B), col = 'black')
legend("topright",inset = 0.05, legend=c("SM1-30A","SM1-30B","SM1-60A","SM1-60B"), col=c("red","orange","blue","black"), lty=1)}

{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 2")
lines(smws32_2mean$Group.1,smws32_2mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM2-30A","SM2-30B","SM2-60A","SM2-60B"), col=c("red","orange","blue","black"), lty=1)}

##extra filtering for smws32_60a
smws32_2mean <- smws32_2mean %>%
  mutate(smois60A = replace(smois60A,Date >= "2017-11-26", NA ))

{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 2 Filtered")
lines(smws32_2mean$Group.1,smws32_2mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM2-30A","SM2-30B","SM2-60A","SM2-60B"), col=c("red","orange","blue","black"), lty=1)}

{plot(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 3")
lines(smws32_3mean$Group.1,smws32_3mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM3-30A","SM3-30B","SM3-60A","SM3-60B"), col=c("red","orange","blue","black"), lty=1)}

soilmoisture30a<- data.frame(smws32_1mean$Date,smws32_1mean$smois30A, smws32_1mean$smois30B,smws32_1mean$smois60A,smws32_1mean$smois60B)
soilmoisture30b<- data.frame(smws32_2mean$Date,smws32_2mean$smois30A, smws32_2mean$smois30B,smws32_2mean$smois60A,smws32_2mean$smois60B)
soilmoisture30c<- data.frame(smws32_3mean$Date,smws32_3mean$smois30A,smws32_3mean$smois30B,smws32_3mean$smois60A,smws32_3mean$smois60B)


merged30a<- merge(soilmoisture30a,soilmoisture30b, by.x = "smws32_1mean.Date", by.y = "smws32_2mean.Date", all = FALSE)

merged30b<- merge(merged30a,soilmoisture30c, by.x = "smws32_1mean.Date", by.y = "smws32_3mean.Date", all = TRUE)

merged30c<- as.data.frame(rowMeans((merged30b[-1]), na.rm= TRUE))


mergedsm<- cbind(merged30b,merged30c)

colnames(mergedsm)<- c("Date","smois30site1a","smois30site1b","smois60site1a","smois60site1b","smois30site2a","smois30site2b","smois60site2a","smois60site2b","smois30site3a","smois30site3b","smois60site3a","smois60site3b","mergedsoilmoisture")

plot(mergedsm$Date, mergedsm$mergedsoilmoisture, type = "l", main = "Average WS32 Soil Moisture", xlab = "Date", ylab = "Merged Soil Moisture All Data")

plot(mergedsm)

## Calculate Calibration/Validation Split based on SM data

daysofsmdata<- as.numeric(range(mergedsm$Date)[2] - range(mergedsm$Date)[1])

#70/30 split

Caldays<- round(0.7*daysofsmdata)
Valdays<- daysofsmdata-Caldays

Caldates<- as.Date(1)

Caldates[1]<- range(mergedsm$Date)[1]

Caldates[2]<- range(mergedsm$Date)[1]+Caldays

Valdates<- as.Date(1)

Valdates[1]<-  range(Caldates)[2]+1

Valdates[2]<- range(mergedsm$Date)[2]

Plot and Filter Observed Streamflow Data

Original Datasets are recorded as Water Year, We convert to calendar year

## import daily flow dataset
obsws32<- read.csv(paste0(filedir,"input_files/streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32.csv"), header = F) 

# dataset is recorded as water year, here we convert from water year and note that the water year for this landscape begins in November
## calendar year seems to be correct

obsws32$WYR <- as.numeric(obsws32$V1)
obsws32$Month <- as.numeric(obsws32$V2)
obsws32$Day <- as.numeric(obsws32$V3)
obsws32$flow <- as.numeric(obsws32$V4)

obsws32$CalendarYear<- ifelse((obsws32$Month=="12"|obsws32$Month=="11"),obsws32$WYR-1,obsws32$WYR)

obsws32$CalendarDate <- as.Date(with(obsws32,paste(CalendarYear,Month,Day,sep="-")),"%Y-%m-%d") 

##this removal of 0s may be reason for lag, 0s may be purposefully included in dataset
# instead create timeseries of date from start to end of timeseries and then merge
#obsws32<- obsws32[!(obsws32$flow=="0"),]

obsws32 <- drop_na(obsws32)

plot(obsws32$CalendarDate, obsws32$flow, type = "l", ylab = "Observed Discharge mm/d", xlab = "Date", main = "Observed Discharge WS32 Older Dataset")

obsws32_1<- read.csv(paste0(filedir,"input_files/streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2017.csv"), header = T)
obsws32_2<- read.csv(paste0(filedir,"input_files/streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2018.csv"), header = T)
obsws32_3<- read.csv(paste0(filedir,"input_files/streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2019.csv"), header = T)
obsws32_4<- read.csv(paste0(filedir,"input_files/streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2020.csv"), header = T)
obsws32_5<- read.csv(paste0(filedir,"input_files/streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2021.csv"), header = T)

obsws32b<- rbind(obsws32_1,obsws32_2,obsws32_3,obsws32_4,obsws32_5)

obsws32b$CalendarDate <- as.Date(with(obsws32b,paste(year,month,day,sep="-")),"%Y-%m-%d")

plot(obsws32b$CalendarDate,obsws32b$flow, type = "l", xlab = "Date", ylab = "Streamflow", main = "Observed Discharge WS32 Newer Datasets Merged")

dataset1   <- data.frame(Date=obsws32$CalendarDate,flow=obsws32$flow)
dataset2   <- data.frame(Date=obsws32b$CalendarDate,flow=obsws32b$flow)

fullset<- rbind(dataset1,dataset2)

### Streamflow at Coweeta logged as CSM
### conversion from EcohydRology project
## because we are converting from CSM we set watershed area to 1

fullset$discharge_mm<- ConvertFlowUnits(cfs = fullset$flow, WA=1, AREAunits = "mi2")

#{plot(fullset$Date,fullset$flow, type = "l", xlab = "Date", ylab = "Flow", main = "Streamflow WS32")
#lines(fullset$Date,fullset$discharge_mm, col = 'red')}

plot(fullset$Date,fullset$discharge_mm, type = "l", xlab = "Date", ylab = "Discharge mm/d", main = "Observed Discharge WS32 All Datasets Merged")

plot(fullset$Date,fullset$discharge_mm, type = "l", xlab = "Date", ylab = " Log Discharge mm/d", main = "Log of Observed Discharge WS32 All Datasets Merged", log = "y")

fullset$observedstreamflow<- fullset$discharge_mm

write.csv(fullset,paste0(filedir,"input_files/streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/ObservedWS32.csv"))

obsws32<- fullset

obsws32sm<- merge(obsws32,mergedsm, by = "Date")

#Caldates
#Valdates

obsws32cal <- obsws32[obsws32$Date >= Caldates[1] & obsws32$Date <= Caldates[2], ]
  
obsws32val <- obsws32[obsws32$Date >= Valdates[1] & obsws32$Date <= Valdates[2], ]
  
obsws32smcal <- obsws32sm[obsws32sm$Date >= Caldates[1] & obsws32sm$Date <= Caldates[2], ]

obsws32smval <- obsws32sm[obsws32sm$Date >= Valdates[1] & obsws32sm$Date <= Valdates[2], ]

Merge Observed and Simulated Data

setwd(system.file("extdata/", package = "RHESSysIOinR"))
getwd()
## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"
singleruntest<- readin_rhessys_output("out/cwws32_runsingle")
## [1] "basin daily and basin daily grow have different number rows"
## merge values instead of subsetting to match up with simulated dates
singleruntest$bd<- merge(singleruntest$bd,obsws32, by.x = "date", by.y = "Date", all = FALSE)

## merge soil moisture values but include all data
singleruntest$data <-merge(singleruntest$bd,obsws32sm, by.x = "date", by.y = "Date", all = TRUE)

Plot Observed and Simulated Data, Plot RHESSys Outputs

## Plot Hydrograph comparing modeled and observed streamflow
{hydrograph(input=singleruntest$bd, streamflow=singleruntest$bd$observedstreamflow, streamflow2 = singleruntest$bd$streamflow, timeSeries = singleruntest$bd$date, precip = singleruntest$bd$streamflow,
           P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')

legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}

{plot(singleruntest$data$date,singleruntest$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Streamflow")
lines(singleruntest$data$date, singleruntest$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}

{plot(singleruntest$data$date,singleruntest$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Log Streamflow", log = "y")
lines(singleruntest$data$date, singleruntest$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}

predictedvsobservedstreamlm<- lm(observedstreamflow~streamflow, data = singleruntest$bd)

{plot(singleruntest$bd$observedstreamflow,singleruntest$bd$streamflow, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Streamflow", xlim = c(0,55), ylim = c(0,55))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedstreamlm)$r.squared,digits=3)))}

{plot(singleruntest$data$date,(singleruntest$data$rz_storage/singleruntest$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(singleruntest$data$date, singleruntest$data$mergedsoilmoisture, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("black","red"), lty=1:1)}

predictedvsobservedsoillm<- lm((rz_storage/rootdepth)~mergedsoilmoisture, data = singleruntest$data)

{plot((singleruntest$data$rz_storage/singleruntest$data$rootdepth),singleruntest$data$mergedsoilmoisture, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Basin Scale Soil Moisture", xlim = c(0.05,0.4), ylim = c(0.05,0.4))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedsoillm)$r.squared,digits=3)))}

## Plot other model outputs

par(mfrow=c(5,1), mar = c(2,5,4,2))

plot(singleruntest$bd$date,singleruntest$bd$tmin, type = "l", col = 'BLUE', ylim = c(-15,30), ylab = "Temperature", xlab = "Date")
lines(singleruntest$bd$date,singleruntest$bd$tmax, type = "l", col = 'RED')
plot(singleruntest$bd$date,singleruntest$bd$precip, type = "h", ylab = "Precipitation", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$streamflow, type = "l", ylab = "Streamflow", xlab = "Date", col = 'blue')
plot(singleruntest$bd$date,singleruntest$bd$baseflow, type = "l", ylab = "Baseflow", xlab = "Date", col = 'brown')
plot(singleruntest$bd$date,singleruntest$bd$rz_storage, type = "l", ylab = "RZ Storage", xlab = "Date", col = 'dark green')

plot(singleruntest$bd$date,singleruntest$bd$psn ,type = "l", ylab = "PSN", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$et, type = "l", ylab = "Evaoptranspiration", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$pet, type = "l", ylab = "PET", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$gw.Qout, type = "l", ylab = "Q out", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$gw.storage, type = "l", ylab = "GW Storage", xlab = "Date")

par(mfrow=c(1,1))

### will only plot in growth mode , used to check that vegetation is initialized properly ##
plot(singleruntest$bd$date, singleruntest$bd$lai, type = "l", main = "LAI", xlab = "Date")

plot(singleruntest$bd$date, singleruntest$bd$soilc, type = "l", main = "Soil Carbon", xlab = "Date")

## plot cal/val sets

{plot(singleruntest$data$date,(singleruntest$data$rz_storage/singleruntest$data$rootdepth), type = "l", ylim = c(0,0.35),xlim = c(obsws32smcal$Date[1],obsws32smval$Date[length(obsws32smval$Date)]), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(obsws32smcal$Date,obsws32smcal$mergedsoilmoisture, col = "BLUE")
lines(obsws32smval$Date,obsws32smval$mergedsoilmoisture, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Soil Moisture","Observed Calibration Soil Moisture", "Observed Validation Soil Moisture"), col=c("black","blue","red"), lty=c(1,1,1))}

{plot(singleruntest$data$date,singleruntest$data$streamflow, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Basin Scale Streamflow",xlim = c(obsws32smcal$Date[1],obsws32smval$Date[length(obsws32smval$Date)]))
lines(obsws32cal$Date,obsws32cal$observedstreamflow, col = "BLUE")
lines(obsws32val$Date,obsws32val$observedstreamflow, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Streamflow","Observed Calibration Streamflow", "Observed Validation Streamflow"), col=c("black","blue","red"), lty=c(1,1,1))}

Testing Precip input vs model precip

Spatial Data Display

Original input maps have been cropped before processing in R, cropping seemed to reduce errors in preparation, below are the input maps used by R

setwd(system.file("extdata/", package = "RHESSysIOinR"))

# Original input maps have been cropped before processing in R, without cropping the RHESSYs preparation in the R environment will not function properly, below are the input maps used by R
ws32<- raster("grassexport/cwt_ws32/basin_ws32.tif")
patches<- raster("grassexport/cwt_ws32/patch.tif")
acc10<- raster("grassexport/cwt_ws32/acc10.tif")
aspect10<- raster("grassexport/cwt_ws32/aspect10.tif")
#bt1000<- raster("grassexport/cwt_ws32/b.t1000.tif")
dem10<- raster("grassexport/cwt_ws32/dem10.tif")
dem10f<- raster("grassexport/cwt_ws32/dem10f.tif")
dir10<- raster("grassexport/cwt_ws32/dir10.tif")
drain10<- raster("grassexport/cwt_ws32/drain10.tif")
#ht1000<- raster("grassexport/cwt_ws32/h.t1000.tif")
hillslope<- raster("grassexport/cwt_ws32/hillslope.tif")
roads<- brick("grassexport/cwt_ws32/roads.tif")
#shade10<- raster("grassexport/cwt_ws32/shade10.tif")
slope10<- raster("grassexport/cwt_ws32/slope10.tif")
streams<- raster("grassexport/cwt_ws32/streams.tif")
topidx10<- raster("grassexport/cwt_ws32/topidx10_100.tif")
#xmap<- raster("grassexport/cwt_ws32/xmap.tif")
#ymap<- raster("grassexport/cwt_ws32/ymap.tif")

soilstatic<- raster("grassexport/cwt_ws32/ssurgosoil.tif")

reclass_df <- c(0,256,9)
reclass_static <- matrix(reclass_df, ncol = 3, byrow = TRUE)
soilstaticmap<- reclassify(soilstatic, reclass_static)


soilssurgo<- raster("grassexport/cwt_ws32/ssurgosoil.tif")
soildsm<- raster("grassexport/cwt_ws32/dsmsoil.tif")
veg<- raster("grassexport/cwt_ws32/nlcdveg.tif")


plot(ws32, main ="ws32")

plot(patches, main = "patches")

plot(acc10, main = "acc10")

plot(aspect10, main = "aspect10")

#levelplot(bt1000, col.regions = brewer.pal(name='Spectral', n = 11),at=seq(0,6500,1))
#plot(bt1000, main = "bt1000")
plot(dem10, main = "dem10")

plot(dem10f, main = "dem10f")

plot(dir10, main = "dir10")

plot(drain10, main = "drain10")

#plot(ht1000, main = "ht1000")
plot(hillslope, main = "hillslope")

plot(roads, main = "roads")

#plot(shade10, main = "shade10")
plot(slope10, main = "slope10")

plot(streams, main = "streams")

plot(topidx10, main = "topidx10")

#plot(xmap, main = "xmap")
#plot(ymap, main = "ymap")

plot(soilstaticmap, main = "Static Soil Input", zlim = c(0,26))

plot(soilssurgo, main = "SSURGO Soil Input", zlim = c(0,26))

plot(soildsm, main = "RSS Soil Input", zlim = c(0,26))

plot(veg, main = "NLCD Veg Input")

New Soil Map Figure

## plot 3 rasters with common legend, legend should only be integers shown in rasters

##IN STEP 1

setwd(system.file("extdata/", package = "RHESSysIOinR"))

# use watershed polygons for BBOX to request elevation data
x <- read_sf(watershedshapefile)
# buffer to ensure that there are no truncated watersheds
x.buff <- st_as_sf(st_buffer(st_as_sfc(st_bbox(x)), dist = 1000))
# convert to GCS WGS84 -> no transformation / warp will be applied to DEM
x.gcs <- st_transform(x, 4326)
x.buff.gcs <- st_transform(x.buff, 4326)
# requires sf collection (geometry + attributes)
# get DEM in GCS WGS84 and use z = 14 for best available data
elevationraster <- get_elev_raster(locations = x.buff.gcs, z = 14, clip = 'bbox')

res(elevationraster)
## [1] 3.921794e-05 3.921794e-05
# check: OK
{plot(elevationraster, main = "DEM, Watershed Outline and Buffer" )
  plot(st_geometry(x.buff.gcs), add = TRUE, col = NA, lwd = 2)
  plot(st_geometry(x.gcs), add = TRUE, col = NA)}

{plot(elevationraster, main = "DEM and selected Watershed")
  plot(st_geometry(x.gcs[which(x$WS=='32'),]), add = TRUE, col = NA)}

cropped_dem<- crop(elevationraster,st_bbox(x.gcs))

plot(cropped_dem, main = "DEM cropped to Buffer")
plot(st_geometry(x.gcs[which(x$WS=='32'),]), add = TRUE, col = NA)

newproj<- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"

reproje<- projectRaster(from = cropped_dem, crs = newproj, method = 'bilinear')

plot(reproje, main = "Reprojected cropped DEM")

dummydata<- reproje

res(dummydata) <- 10

resampleddem <- resample(reproje, dummydata, method = 'bilinear')

plot(resampleddem, main = "DEM reprojected to 10m Resolution")

{plot(resampleddem, main = "DEM reprojected to 10m Resolution")
  plot(st_geometry(x.buff.gcs), add = TRUE, col = NA, lwd = 2)
  plot(st_geometry(x.gcs), add = TRUE, col = NA)}

##In STEP1 Quick map creation for WS7 soil moisture data display

ws7demcrop <- crop(resampleddem,x[27,])
ws7demcrop <- mask(ws7demcrop,x[27,])
plot(ws7demcrop)

Observedws7SMLocationsnt <- read_sf("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1305.19/1305.kml")
Observedws7SMLocations <- st_transform(Observedws7SMLocationsnt, crs(dem10))

Newsoilmoisturelocationsnt<- read_sf(newsmlocations)
Newsoilmoisturelocations<- st_transform(Newsoilmoisturelocationsnt, crs(dem10))

colorramp<- viridis(1028)

{plot(ws7demcrop, col = colorramp, main = "Coweeta Observed Soil Moisture Locations WS32", legend.args=list(text="Elevation(m)", side = 3, line = 2))
points(c(Observedws7SMLocations[[3]][[3]][1],Observedws7SMLocations[[3]][[4]][1],Observedws7SMLocations[[3]][[5]][1]), c(Observedws7SMLocations[[3]][[3]][2],Observedws7SMLocations[[3]][[4]][2],Observedws7SMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)
points(st_coordinates(Newsoilmoisturelocations), pch =21, col = "black", bg = "green", lwd = 1)
legend("topright", legend = c("Legacy","New"), pch = 20, col= c("red","green"), cex = 1.2)
}

Show Observed Soil Moisture Locations

ObservedSMLocationsnt <- read_sf(paste0(filedir,"input_files/soil_moisture/knb-lter-cwt.1308.19/1308.kml"))

ObservedSMLocations <- st_transform(ObservedSMLocationsnt, crs(dem10))

Newsoilmoisturelocationsnt<- read_sf(newsmlocations)

Newsoilmoisturelocations<- st_transform(Newsoilmoisturelocationsnt, crs(dem10))

colorramp<- viridis(1028)

{plot(dem10, col = colorramp, main = "Coweeta Observed Soil Moisture Locations")
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)

}

 {plot(dem10, col = colorramp, main = "Coweeta Observed Soil Moisture Locations WS32", legend.args=list(text="Elevation(m)", side = 3, line = 2))
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)
points(st_coordinates(Newsoilmoisturelocations), pch =21, col = "black", bg = "green", lwd = 1)
legend("topright", legend = c("Legacy","New"), pch = 20, col= c("red","green"), cex = 1.2)
}

ObservedSMLocations$Name[3]
## [1] "132 plot center"
ObservedSMLocations[[3]][[3]]

ObservedSMLocations$Name[4]
## [1] "232 plot center"
ObservedSMLocations[[3]][[4]]

ObservedSMLocations$Name[5]
## [1] "332 plot center"
ObservedSMLocations[[3]][[5]]

matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 3)
##      [,1]              [,2]               [,3]              
## [1,] "132 plot center" "275432.649744231" "3881522.05364127"
## [2,] "232 plot center" "275578.214010129" "3881428.55011708"
## [3,] "332 plot center" "275760.864807773" "3881284.18476691"
library(mapview)

Coweeta <- st_read(watershedshapefile)
## Reading layer `Coweeta_Hydrologic_Laboratory' from data source 
##   `C:\Users\Carlos\Documents\GitHub\RHESSys74-R-Markdown\input_files\shapefiles\watersheds\Coweeta_Hydrologic_Laboratory.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 273971.1 ymin: 3877997 xmax: 279140.9 ymax: 3883954
## Projected CRS: NAD83 / UTM zone 17N
mapview(Coweeta, label = Coweeta$WS, alpha.regions = 0.1, alpha = 1) +
  mapview(ObservedSMLocations[3:5,]) 
par(mfrow=c(1,2),mar=c(2, 3, 2, 5))
 {plot(dem10, col = colorramp, main = "WS32", legend.args=list(text="          Elevation(m)", side = 3, line = 2, cex = 0.8))
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)
points(st_coordinates(Newsoilmoisturelocations), pch =21, col = "black", bg = "green", lwd = 1)
legend("topright", legend = c("Legacy","New"), pch = 20, col= c("red","green"), cex = 1.2)
}

{plot(ws7demcrop, col = colorramp, main = "WS7", legend.args=list(text="          Elevation(m)", side = 3, line = 2, cex = 0.8))
points(c(Observedws7SMLocations[[3]][[3]][1],Observedws7SMLocations[[3]][[4]][1],Observedws7SMLocations[[3]][[5]][1]), c(Observedws7SMLocations[[3]][[3]][2],Observedws7SMLocations[[3]][[4]][2],Observedws7SMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)
points(st_coordinates(Newsoilmoisturelocations), pch =21, col = "black", bg = "green", lwd = 1)
legend("topright", legend = c("Legacy","New"), pch = 20, col= c("red","green"), cex = 1.2)
}

par(mfrow=c(1,1))

Display Observed and Simulated Soil Moisture

## Extract values for single date
displayday <- singleruntest$pd[which(singleruntest$pd$date=="2018-10-21")]

## reclassify all patches to na and then reclassify to model output values
reclass_all<- c(-2147483648, 2147483647,NA)
reclass_allm<- matrix(reclass_all, ncol = 3, byrow = TRUE)
displaydaypatches<- reclassify(patches,reclass_allm)

displayday$calculatedrzsm <- displayday$rz_storage/displayday$root.depth

reclass_displayday <- cbind(displayday$patchID,displayday$calculatedrzsm)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches<- reclassify(patches,reclass_displayday)

## Plot Spatial Data Output for RZ Storage
#plot(displaydaypatches, xlim = c(280900,282500), ylim = c(4870000,4872000))

{plot(mask(displaydaypatches,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 RZ Soil Moisture Prediction for 1 day", line = -2)
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", cex = 2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$root.depth)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches2<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches2,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 root.depth Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$rz_storage)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches3<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches3,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 rz_storage Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$sat_def)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches4<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches4,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 saturation deficit Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$rz_field_capacity/1000)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches5<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches5,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 RZ Field Capacity for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

## Plot Hydrograph
{hydrograph(input=singleruntest$bd, streamflow=singleruntest$bd$observedstreamflow, streamflow2 = singleruntest$bd$streamflow, timeSeries = singleruntest$bd$date, precip = singleruntest$bd$streamflow,
           P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')
legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}

### end show patch map

patchnumbers<- raster::extract(patches,matrix(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 2))

ObserveSoilMoisturePatches <- matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2],patchnumbers), nrow = 3, ncol = 4)

## Extract values for a single patch, select patch ID
ObserveSoilMoisturePatches[1,4]
## [1] "136557"
displaypatch<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[1,4])]
ObserveSoilMoisturePatches[2,4]
## [1] "141942"
displaypatch2<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[2,4])]
ObserveSoilMoisturePatches[3,4]
## [1] "149478"
displaypatch3<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[3,4])]

plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #147322', ylim = c(0,.4))

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", ylim = c(0,.4), col = 'darkorange', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patches",
     xlab = "Date", ylab = "Soil Moisture mm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A,type = "l", col = 'darkred')
lines(smws32_3mean$Group.1,smws32_3mean$smois30A,type = "l", col = 'blue')
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'red', lty = 2)
lines(displaypatch2$date,(displaypatch2$rz_storage/displaypatch2$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'green', lty = 2)
lines(displaypatch3$date,(displaypatch3$rz_storage/displaypatch3$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'black', lty = 2)

legend("bottomleft",inset = .01, legend=c("Observed Patch 1","Predicted Patch 1","Observed Patch 2","Predicted Patch 2","Observed Patch 3","Predicted Patch 3"), col=c("darkorange", "red","darkred","green","blue","black"), lty=c(1,2,1,2,1,2))}

{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", ylim = c(0,0.4), col = 'BLUE', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patch #154914",
     xlab = "Date", ylab = "Soil Moisture mm")
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #154914', col = 'red')
legend("bottomleft",inset = .01, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("blue","red"), lty=1:1)}

#obsws32cal <- obsws32[obsws32$Date >= Caldates[1] & obsws32$Date <= Caldates[2], ]
#obsws32val <- obsws32[obsws32$Date >= Valdates[1] & obsws32$Date <= Valdates[2], ]
#obsws32smcal <- obsws32sm[obsws32sm$Date >= Caldates[1] & obsws32sm$Date <= Caldates[2], ]
#obsws32smval <- obsws32sm[obsws32sm$Date >= Valdates[1] & obsws32sm$Date <= Valdates[2], ]

head(obsws32cal)
##             Date flow discharge_mm observedstreamflow
## 22961 2017-01-01 2.30     2.172645           2.172645
## 22962 2017-01-02 3.04     2.871669           2.871669
## 22963 2017-01-03 2.72     2.569388           2.569388
## 22964 2017-01-04 2.11     1.993165           1.993165
## 22965 2017-01-05 1.92     1.813686           1.813686
## 22966 2017-01-06 1.81     1.709777           1.709777
head(obsws32val)
##             Date flow discharge_mm observedstreamflow
## 23311 2017-12-17 2.14     2.021504           2.021504
## 23312 2017-12-18 2.27     2.144306           2.144306
## 23313 2017-12-19 2.23     2.106521           2.106521
## 23314 2017-12-20 3.39     3.202289           3.202289
## 23315 2017-12-21 2.74     2.588281           2.588281
## 23316 2017-12-22 2.60     2.456033           2.456033
head(obsws32smcal)
##         Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 1 2017-01-01 2.30     2.172645           2.172645     0.2194375     0.2468125
## 2 2017-01-02 3.04     2.871669           2.871669     0.2241528     0.2509896
## 3 2017-01-03 2.72     2.569388           2.569388     0.2184132     0.2409236
## 4 2017-01-04 2.11     1.993165           1.993165     0.1952743     0.2236250
## 5 2017-01-05 1.92     1.813686           1.813686     0.1825451     0.2149410
## 6 2017-01-06 1.81     1.709777           1.709777     0.1745729     0.2097500
##   smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 1     0.1988021     0.2463715     0.2162465     0.2141181     0.2191528
## 2     0.2243646     0.2737917     0.2244410     0.2392847     0.2239028
## 3     0.2177813     0.2657361     0.2209722     0.2291771     0.2154028
## 4     0.1993090     0.2483194     0.2058437     0.2039549     0.1967431
## 5     0.1845556     0.2335799     0.1972882     0.1883021     0.1870486
## 6     0.1744132     0.2231181     0.1914861     0.1774236     0.1808646
##   smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 1     0.1775868     0.2830729     0.2573576     0.2625312     0.2399167
## 2     0.2186528     0.2833924     0.2626007     0.2942917     0.2684410
## 3     0.2180486     0.2697604     0.2480833     0.2926944     0.2704132
## 4     0.1968403     0.2399167     0.2206285     0.2794965     0.2609514
## 5     0.1786215     0.2264861     0.2081111     0.2668299     0.2512847
## 6     0.1669340     0.2187014     0.2009896     0.2569028     0.2437292
##   mergedsoilmoisture
## 1          0.2317839
## 2          0.2490255
## 3          0.2422839
## 4          0.2225752
## 5          0.2099661
## 6          0.2015738
head(obsws32smval)
##           Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 351 2017-12-17 2.14     2.021504           2.021504     0.1181493     0.2395313
## 352 2017-12-18 2.27     2.144306           2.144306     0.1165625     0.2392708
## 353 2017-12-19 2.23     2.106521           2.106521     0.1135139     0.2332326
## 354 2017-12-20 3.39     3.202289           3.202289     0.1484514     0.2609861
## 355 2017-12-21 2.74     2.588281           2.588281     0.1391528     0.2532083
## 356 2017-12-22 2.60     2.456033           2.456033     0.1312049     0.2428194
##     smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 351     0.1677882     0.2054826     0.1575382     0.1352778            NA
## 352     0.1635625     0.2064132     0.1605903     0.1370000            NA
## 353     0.1600174     0.2057812     0.1585660     0.1363750            NA
## 354     0.1844479     0.2319132     0.1872778     0.1721458            NA
## 355     0.2092639     0.2517917     0.1925139     0.1736111            NA
## 356     0.1954965     0.2400937     0.1847604     0.1649479            NA
##     smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 351    0.09663194     0.2347188     0.1869757     0.2350000     0.2169896
## 352    0.09784375     0.2385694     0.1832847     0.2350000     0.2167847
## 353    0.09828472     0.2289757     0.1792153     0.2350000     0.2160000
## 354    0.12405208     0.2609826     0.2203681     0.2421354     0.2173993
## 355    0.14059028     0.2525278     0.2136979     0.2712326     0.2365937
## 356    0.13011806     0.2397049     0.1995382     0.2683160     0.2404201
##     mergedsoilmoisture
## 351          0.1812803
## 352          0.1813529
## 353          0.1786329
## 354          0.2045600
## 355          0.2121985
## 356          0.2034018
#subset Validation range based on Val Dates set by available SM data

validationsubset <- singleruntest$bd[singleruntest$bd$date >= Valdates[1] & singleruntest$bd$date <= Valdates[2], ]

## Fit tests for streamflow

NSE(validationsubset$streamflow,validationsubset$observedstreamflow)
## [1] 0.2504567
NSE(validationsubset$streamflow,validationsubset$observedstreamflow, FUN=log, epsilon = "Pushpalatha2012", na.rm=TRUE)
## [1] 0.2194174
nse_value <- NSE(log(validationsubset$streamflow + 1), log(validationsubset$observedstreamflow + 1), na.rm = TRUE)


 KGE(validationsubset$streamflow,validationsubset$observedstreamflow)
## [1] 0.4080325
plot(validationsubset$streamflow~validationsubset$date, type = "l", col = 'red')
lines(validationsubset$observedstreamflow~validationsubset$date)

# Plot Observed and Predicted Streamflow

{plot(validationsubset$date,validationsubset$streamflow, type = 'l', lty = 2, col = 'RED', main = "Streamflow")
lines(validationsubset$date,validationsubset$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}

{plot(singleruntest$bd$date,singleruntest$bd$streamflow, type = 'l', lty = 2, col = 'RED', main = "Streamflow")
lines(singleruntest$bd$date,singleruntest$bd$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}

#Plot Model RZ SM outputs for patch
{plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l",ylim = c(0,0.50), ylab = "Root Zone Storage (mm/m) ", xlab = "Date", main = 'Root Zone Display Patch')
lines(displaypatch$date,(displaypatch$rz_field_capacity/displaypatch$root.depth), type = "l", col = "blue")
lines(displaypatch$date,(displaypatch$rz_wilting_point/displaypatch$root.depth), type = "l", col = "red")
legend("bottomleft",inset = .1, legend=c("Predicted RZ Storage","RZ Field Capacity","RZ Wilting Point"), col=c("black","blue","red"), lty=1:1)}